Graphene buffer layer on Si-terminated SiC : a mult i- minima energy surface studied 

with an empirical interatomic potential 



Evelyne Lampin,* Catherine Priester, and Christophe Krzeminski 
IEMN - BP 60069 - 59652 Villeneuve d'Ascq Cedex - France 

Laurence Magaud 
Institut Neel - BP 166 - 38042 Grenoble Cedex 9 - France 
(Dated: December 10, 2009) 

The atomistic structure of the graphene buffer layer on Si-terminated SiC is studied using a modi- 
fied environment-dependent interatomic potential (EDIP). The investigation of equilibrium state by 
conjuguate gradients suffers from a complex multi-minima energy surface. A dedicated procedure 
is therefore presented to provide a suitable initial configuration on the way to the minimum. The 
result forms an hexagonal pattern with unsticked rods to release the misfit with the surface. The 
structure presents an agreement with the global pattern obtained by experiments and even with the 
details of an ah initio calculation. 

I. INTRODUCTION 

Graphene is an attractive alternative to silicon for high- mobility devices due to its unique electronic properties. 
The graphitisation of SiC by high-temperature annealing 1 is often seen as the most promising route to fabricate large 
layers of graphene. Although many experimental works were devoted to the refinement of the fabrication process, 
the graphitization mechanism remains poorly understood. In particular, the control of the sublimation of Si atoms 
during and after formation of the first graphene layer is a key step of the fabrication of larger and more regular 
layers 2 , ab initio calculations have shown that the first graphene layer, the so-called buffer layer, on the Si-terminated 
surface has not a graphenelike dispersion 3 and is strongly bonded to the substrate 4 , 5 . The experimental counterpart 
is not simple, in particular due to the high sensitivity to the preparation procedure 6 . Further understanding could be 
provided with molecular dynamics calculations that describe the evolution of an atomic system submitted to thermal 
annealing. Such calculations could only be carried out with rough approximations in the framework of ab initio 
descriptions and the use of an interatomic potential is therefore more adapted. We present in this study the results of 
calculations performed with a version of the environment-dependent interatomic potential (EDIP) modified for SiC 
by Lucas 7 . The equilibrium state of the buffer layer on SiC (0001) is studied. The difficulty to achieve convergence 
in this multi-minima energy surface is demonstrated, and a specific procedure is given to reach equilibrium. The 
configuration of minimal energy presents the symetries observed by experiments. Moreover, the details are compared 
to ab initio results. A good agreement is obtained and the optimum even presents an higher level of order. We 
therefore conclude on the ability of the potential to describe the graphitisation of SiC surface. 

II. METHOD 

In a first step, a model of graphene on SiC(0001) is constructed. For the sake of simplicity we consider a 3C-SiC 
substrate, but the system is perfectly equivalent to 4H- and 6H-SiC provided a planar SiC surface is considered (the 
equivalence does not hold in the case of a step for example). The lattices of graphene and SiC do not match. A low 
misfit is obtained for a (6\/3 x 6\/ r 3)R30° supercell 8 put in regard of a sheet of 13-hexagon wide honeycomb. The 
misfit, defined by the ratio : 

13 X ttgraphene 

6^ x a Si c 

is equal to 0.7 % experimentally and 3.0 % with the lattice parameters (a gra phene 7 &Sic) obtained at equilibrium using 
the interatomic potential presented hereafter. The graphene layer of the basic cell contains 338 atoms of C. The depth 
of the SiC substrate has to be high enough to absorb the strain generated by misfit. We choose to use 6 SiC bilayers, 
with a total of 1296 atoms, although the dependence of the results on the depth of the substrate will be tested in the 
last section of this article. The basic hexagonal cell is replicated twice along x and y to facilitate the visualization of 
the patterns formed by atoms in the buffer layer. Fig. 1 presents lateral and top views of the corresponding sytems 
of 6536 atoms. 

The interactions between atoms are described by an empirical interatomic potential. Empirical interatomic po- 
tentials are developed for a specific domain of a material physics. It is generally difficult to estimate the range of 
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FIG. 1: Lateral (left) and top (right) views of the atomic system under study. The C atoms are colored in orange, the Si in 
grey. The delimitations of the basic cell are shown on the right view. 



validity of a given potential when it is applied to a different phenomenon for the first time 9 . Indeed one of the aim of 
this work is to confront the potential under study to ab initio calculations in order to validate its application to the 
study of graphene on SiC. The interatomic potential is a slightly modified version of EDIP 10 , 11 , 12 parametrised for 
C-C and Si-C interactions by Lucas 7 . Other classical interatomic potentials such as Tersoff 13 and Brenner 14 have also 
been adapted and parametrised for Si-C 15 , 16 , 17 , 18 , 19 . Tewary and Yang 20 have constructed a parametric interatomic 
potential but it is restricted to graphene (only C-C interations). As for the study of SiC graphitization, Tang et al. 
used Tersoff 's potential in molecular dynamics simulation of the heating of Si-less layers on SiC 21 and of graphene 
nanoribbons 22 although the coordination of order 3 is not natural for a potential developed for diamond lattice. In the 
present study, the environment-dependent mathematical formulation of EDIP has been assumed to be an advantage 
for the modelisation of the bonding of graphene layer on the SiC surface including sp2 and sp3 bonds. The potential 
energy of a system of atoms at position fi is given by : 



i 

Ei = ^2v 2 (r i ,r j ,Z i ) + ^ V$(r u fj, r k , Zi). 



(2) 



The coordination number is : 



exp 



o, 



r < c 
A(r), c < r < a 
r > a. 



(3) 



The corrective function A(r) is added by Lucas 7 to avoid the second nearest neighbour Si-Si interactions in SiC to be 
accounted for. The function is equal to unity except if two Si atoms are interacting : 



exp 



A(r) = < 

with 5 = 0.3 A. The two-body part is given by : 



«/(i-(^) 3 ) 



V 2 (f i ,f j ,Z i ) = A 
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0, 



r < c 
, c < r < a — 5 
r > a — 5 



(4) 
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A(ry)- 
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The three-body term is made of radial and angular contributions : 



Vz{ri,rj,r k ,Zi) = exp ( — - — J A(r ij )exp ( — - — ) A(r ik )h(cos9 ijk , Zi), 



2 



h(l, Z) = X [(1 - exp[-Q(Z)(Z + t(Z)) 2 }) + rjQ(Z)(l + r(Z)) 2 ] , 
Q(Z) = Q exp{-fiZ), t(Z) = u x + u 2 (u 3 e- u * z - e~ 2u * z ), 



(6) 
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Si-Si 


C-C 


Si-C 


a(A) 


3.1213820 


2.2918049 


2.6746212 


c(A) 


2.5609104 


1.5698566 


2.0050591 


a 


3.1083847 


1.1268088 


2.1175967 


A(eV) 


7.9821730 


11.6773355 


9.6545591 


B(A) 


1.5075463 


0.9612016 


1.2037674 


P 


1.2085196 


2.8779528 


2.0432362 


P 


0.0070975 


0.0255960 


0.0163467 




0.5774108 


0.7103453 


0.6404382 


7(A) 


1.1247945 


1.2258251 


1.1742237 



TABLE I: Two-body parameters of the EDIP potential used in the present work 12 , 7 . 





Si-Si-Si 


C-C-C 


Si-SiC 


C-SiC 


Si-CC 








or -CSi 


or -CSi 


or C-SiSi 


A(eV) 


1.4533108 


1.7800184 


1.5349877 


1.6983415 


1.6166646 


rj 


0.2523244 


0.8115872 


0.3921401 


0.6717715 


0.5319558 


Qo 


312.1341346 417.5476826 338.4875216 391.1942956 364.8409086 




0.6966326 


0.6122148 


0.6755281 


0.6333192 


0.6544237 



TABLE II: Three-body parameters of the EDIP potential used in the present work 12 , 7 . 

with u\ = —0.165799, ui = 32.557, us = 0.286198 and 114 = 0.66. The parameter set is given in Tables I and II. The 
result is a potential able to describe Si- and C-based crystals and their defects. The only noticeable drawback, as in 
standard ab initio calculations, is that the Van der Waals interactions are poorly described. Therefore graphite could 
only be obtained by fixing the distance between two planes of graphene. 

Periodic boundary conditions are applied in the lateral directions x and y (see Fig. 1. The bottom SiC bilayer in 
the z direction is frozen to avoid reconstruction of the C-terminated surface. The relaxation is performed with the 
method of the conjugate gradients 23 until a precision of 10 -8 is reached on the total energy (tenths of /ieV/atom) 
after tenths to hundreds of steps. The initial gap between graphene and the surface is set between 1.5 and 3 A by 
step of 0.01 A, the cutoff of the Si-C interactions being equal to 2.67 A and the equilibrium Si-C distance to 1.91 A. 
Three lateral locations of the graphene plane compared to the SiC surface are tested : 

• shift 1 : a C atom in graphene at the vertical of a Si surface atom 

• shift 2 : center of a graphene hexagon at the vertical of a Si surface atom 

• shift 3 : middle of a C-C bond in graphene at the vertical of a Si surface atom. 

III. RESULTS 

The total energy after relaxation is presented in Fig. 2 as a function of the gap between the plane of graphene and 
the surface and for the three shifts. The second and third shifts, corresponding to an alignment on the center of an 
hexagon or on the middle of a bond, give exactly the same energy. Indeed these configurations should be equivalent. 
In the three cases, the result strongly depends on the initial gap. Above the Si-C cutoff radius of the interactions 
(2.67 A), the energy does not depend on the initial gap as expected and corresponds to the sum of the independent 
contributions of the SiC substrate and of the plane of graphene. In the range [2.1-2.6]A, the energy only slighlty 
evolves with the initial gap while there is a strong variation of the result when the graphene plane is set at less than 
2.1 A of the SiC surface. From one value of the gap to the neighboring one, the bonds between graphene and Si 
surface involve different atoms and this is enough to break bonds and form new ones. This is the explanation of the 
great dependence of the total energy on the initial gap of graphene, and the presence of multiminima in energy. The 
method of conjuguate gradient is obviously not able to allow the system to go out of a local minimum. The open 
question is therefore whether the absolute minima is found or not with this method. 

Moreover, we studied in more details the structure of the graphene plane and its bonding to the surface for the 
two gaps that give a low energy in the configuration of shift 1, i. e. an alignment of a C in graphene on a Si at the 
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FIG. 2: Total potential energy versus initial gap between graphene and surface. 



surface. These gaps are 1.63 and 2.26 A. Color maps of the position of the C atoms along z superimposed on dots 
that figure the C position are shown on Figs. 3 and 4 respectively. The patterns have an order 3 symmetry, although 




FIG. 3: Color map of the position of the C atoms in 

graphene along z (black dots). Initial gap set to 1.63 A. FIG. 4: Same figure than Fig. 3 but for an initial gap of 
Blue (red) means near (far from) the SiC surface. Light 2.26 A. 
crosses indicate the position of the Si atoms of the surface. 

an order 6 could be expected, and the structure is surprisingly not regular in Fig. 3. Clearly, the minimisation does 
not seem achieved. The two figures however evidence the complexity of the energy surface, since the topologies of the 
graphene layer are so different with total energies almost equal (difference of 0.007%). 

Therefore we decided to refine the construction of the initial structure to facilitate the relaxation. Indeed, it seems 
obvious that the multiple possible binding of some of the C atoms of the graphene layer to the Si of the surface is 
too difficult to handle by the procedure of conjuguate gradients. Since the choice of the initial configuration is of 
primary importance, 3 initial structures have been specifically built as illustrated in Figs. 5 to 7 In these structures, 
the graphene layer is no longer planar, but the C atoms are set at 3 different heights z. The C atoms initially close to 
a Si at the surface are set at the lower z to strengthen the bonds, their immediate surrounding to an intermediate z 
and the remaining atom to a higher z. By varying the initial mean gap between graphene and the surface, the layer 
is more or less blue, i.e. more or less sticked to the SiC surface. Fig. 5 corresponds to a mean gap of 2.55 A, Fig. 6 
to 2.58 A and Fig. 7 to 2.62 A. A cross view of the system shown on Fig. 5 is also given in Fig. 8 and evidences the 
slight separation along z of the three levels. 

After relaxation, the graphene sheet presents a pronounced buckling of magnitude typically equal to 1.5 A (Fig. 
9). The energy obtained with the three 3- levels systems is smaller than the minimum obtained using a planar sheet 
as a starting point (Fig. 2). The most stable structure is obtained using the highly-sticked system presented in Fig. 
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FIG. 5: Strongly sticked three-level structure of the 

graphene layer. Blue dots are for C atoms closed to the FIG. 6: Moderately sticked three-level structure of the 
surface, red dots for those far from the surface and green graphene layer. Mean gap equal to 2.58 A. Same color 
ones for intermediate values of z. Mean gap equal to 2.55 code than in Fig. 5. 
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FIG. 7: Weakly sticked three-level structure of the 
graphene layer. Mean gap equal to 2.62 A. Same color 
code than in Fig. 5. 
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FIG. 8: Cross view of the system presented in Fig. 5. 



FIG. 9: Cross view after relaxation of the system pre- 
sented in Fig. 5. 



5. A color map of the positions along z of the C atoms in graphene is given in Fig. 10 for the corresponding system. 
Low regions form circular and triangular patterns separated by rods at higher z. These rods are necessary to release 
the misfit of 3.0 % between the graphene sheet and the substrate. The duplication of the elementary pattern results 
in the formation of an hexagonal structure formed of a central disk surrounded by six triangles separated by rods 
distributed in a six-branch star. The cross view shown in Fig. 9 suggests that the C atoms in these rods are not 
bonded to the surface, an assumption that will be discussed in the next section. By comparing the initial structure 
(Fig. 5) to the relaxed one, we guess that the half rods already present initially extend to form complete rods during 
the relaxation. 



IV. DISCUSSION 



The result presented in Fig. 10 has several interesting features. First, it has the lowest energy ever found in our 
various calculations. Moreover, the structure presents a high level of order, in particular compared to solutions such as 
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FIG. 10: Color map of the positions along z of the C atoms (black dots) in graphene. Most stable configuration. Blue (red) 
means near (far from) the SiC surface. Light crosses indicate the position of the Si atoms of the surface. 





ab initio 


Empirical potential 


Interatomic interations 


GGA 24 (code VASP 25 ) 


Modified version of EDIP 7 


Misfit graphene/SiC 


0.7% 


3.0 % 


Basic cell 


1310 atoms, 4 SiC bilayers, 


6536 atoms, 6 SiC bilayers, 


characteristics 


H-passivation at bottom 


frozen atoms at bottom 


Convergence 


Residual forces 


Variation in energy 


criterium 


< 15 meV/A 


< tenths of /xeV/atom 


Number of steps 




Tenths to 


for relaxation 




hundreds 



TABLE III: Main features of the ab initio and empirical potential calculations. 



the one presented in Fig. 3. The symetry level is also the highest, and the (6-\/3 x 6^3)R30° and (6 x 6) periodicities 
already observed by scanning tunneling microscopy 6 are visible in Fig. 10. This solution presents all the features 
expected for the equilibrium state. 

The gap between the SiC surface and the graphene layer, defined in this case as the distance between the upper 
Si atom of the surface and the lower C of the graphene, is equal to 1.68 A. This gap is smaller than the equilibrium 
distance of the Si-C bond with EDIP, i. e. 1.91 A, which is due to the formation of leaning connections between atoms 
that are not at the vertical of each other. The higher C atoms in graphene are at z = 3.17 A, a value to compare to 
the cutoff radius of the Si-C interactions with EDIP equal to 2.67 A. As foreseen in the previous section, the C atoms 
in the rods are indeed not bonded to the surface, their location is only fixed by sp2 interactions with their neighbours 
in the buffer layer. The amplitude of graphene buckling is equal to 1.50 A. The energy per C atom in graphene is 
equal to -7.42 eV. This value is lower than the cohesive energy in graphite (-7.38 eV), which is due to the formations 
of Si-C bonds. Moreover, the impact of the substrate thickness, which could be important due to the stress of the 
bonded graphene in misfit has been studied. The test has been performed on the energy per graphene atom since the 
total number of atom is not constant. It has been obtained a variation of only 0.0056% from 4 to 6 SiC bilayers, and 
0.0005 % from 6 to 15 bilayers. Obviously 4 SiC bilayers are enough to absorb the strain. 

The present result obtained with the EDIP empirical potential is afterwards compared to the result of an ab initio 
calculation 4 (Table III). The structure is duplicated and a color map is applied to visualize the variations in height 
of the C in graphene as in Fig. 10. The result is presented in Fig. 11. The structure is less ordered, and an 
elementary pattern does not clearly stand out, although the hexagonal structure is undoubtly appearent at a larger 
scale. The rods are smaller and the low regions are not so regular but circular forms and quasi triangular ones can 
be distinguished. The difficulty to find the minimum is probably here again at the origin of the irregularities of the 
structure, but the pattern globally corresponds to the result of the EDIP calculation. The atomic positions along z 
and the corresponding Fourier transforms are drawn at larger scale in Fig. 12. The Fourier transforms evidence some 
particular features in common (hexagon at 0.3 A -1 for the 6x6-SiC modulation and 0.7 A -1 for a local symmetry of 
the more or less circular grains at low z) and the higher level of symmetry obtained with the present minimisation 
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FIG. 11: Color map of the positions along z of the C atoms (black dots) in graphene. Result of the ab initio calculation 
presented in Ref. 4. Blue (red) means near (far from) the SiC surface. 

(second hexagon at 0.3 A -1 coming from the 30°-rotated 6x6-SiC periodicity and central hexagon at 0.2 A -1 due to 
the (6a/3 x 6\/3)R30 o elementary pattern, at least more intense in the present calculation). 

The buckling amplitude of the graphene layer is equal to 1.2 A to be compared to 1.5 A for EDIP, for respectively 
0.7 and 3.0 % of graphene-SiC misfit. The space between graphene and the surface is equal to 1.8 A, versus 1.7 A for 



The global agreement of the two calculations leads us to the conclusion that this empirical potential is perfectly 
suited for the study of graphene on SiC. The level of order of the superstructure is even more pronounced, due to the 
methodology developed to find the minimum despite an empirical physical description of atomic interactions. If in 
this case, ab initio calculations were still affordable, the use of an empirical description will be essential to study larger 
systems with surface steps or perform molecular dynamics simulations at a given temperature in order to investigate 
the mechanism of Si sublimation that results in the formation of graphene. 



In conclusion, the bonding of the graphene buffer layer to the Si-terminated surface of SiC has been studied using 
a modified version of the interatomic potential EDIP. The difficulty to find the equilibrium configuration has been 
emphasized, and a dedicated procedure has been developed in order to produce a starting configuration on the 
way to the minimum. The result presents high level of order and symmetry, with an hexagonal pattern formed of 
sticked circular and triangular regions and unbonded rods to release the misfit. The structure is coherent with the 
result obtained using an ab initio description and presents all the symetries identified in experiments. The physical 
description is therefore validated for further study on the graphitization process. 
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FIG. 12: Top: color map on the poisitons of C in graphene. Bottom: Fourier transform of the top color map, zoom at low k 
to emphasize the large scale order. Left: EDIP. Right: ab initio 4 . 
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